void run_sim(Int_t nEvents = 30, TString mcEngine = "TGeant4", Bool_t isMT = false)
{

    // Output file name
    TString outFile = "test.root";

    // Parameter file name
    TString parFile = "params.root";

    // -----   Timer   --------------------------------------------------------
    TStopwatch timer;
		timer.Start();
		// ------------------------------------------------------------------------
		FairLogger *logger = FairLogger::GetLogger();
		// define log file name
		logger->SetLogFileName("MyLog.log");
		// log to screen and to file
		logger->SetLogToScreen(kTRUE);
		logger->SetLogToFile(kTRUE);
		// Print very accurate output. Levels are LOW, MEDIUM and HIGH
		// logger->SetLogVerbosityLevel("HIGH");
		logger->SetLogVerbosityLevel("LOW");
	
    // -----   Create simulation run   ----------------------------------------
    FairRunSim* run = new FairRunSim();
    run->SetName(mcEngine);                        // Transport engine
    run->SetSink(new FairRootFileSink(outFile));   // Output file
    run->SetIsMT(isMT);
    FairRuntimeDb* rtdb = run->GetRuntimeDb();
    // ------------------------------------------------------------------------
		run->SetWriteRunInfoFile(kFALSE);
    // -----   Create media   -------------------------------------------------
    run->SetMaterials("media.geo");   // Materials
    // ------------------------------------------------------------------------

    // -----   Create geometry   ----------------------------------------------

    FairModule* cave = new EtfCave("CAVE");
    cave->SetGeometryFileName("cave.geo");
    run->AddModule(cave);

		FairModule* magnet = new EtfMagnet("Magnet");
		// magnet->SetGeometryFileName("magnet.geo");
		magnet->SetGeometryFileName("magnet_simple.root");
		run->AddModule(magnet);

    FairModule* pipe = new EtfPipe("Pipe");
//    run->AddModule(pipe);

    FairDetector* tofw= new EtfTofw("EtfTofw", kTRUE);
		// tofw->SetGeometryFileName("tofw_simple.root");
		tofw->SetGeometryFileName("tofw_v1.root");
		run->AddModule(tofw);

    FairDetector* neutronWall= new EtfNeutronWall("EtfNeutronWall", kTRUE);
		neutronWall->SetGeometryFileName("neutronWall_simple.root");
		run->AddModule(neutronWall);

    FairDetector* mwdc= new EtfMwdc("EtfMwdc", kTRUE);
		//mwdc->SetGeometryFileName("mwdc_simple.root");
		mwdc->SetGeometryFileName("mwdc_v1.root");
		run->AddModule(mwdc);
    // ------------------------------------------------------------------------

    // -----   Magnetic field   -------------------------------------------
    // Constant Field
		EtfConstField* fMagField = new EtfConstField();
		// FairConstField* fMagField = new FairConstField();
		fMagField->SetFieldRegion(-40,40,-14,14,-51,51);   // values are in cm
    // fMagField->SetFieldRegion(-100,100,-100,100,-100,100);   // values are in cm
    // fMagField->SetField(0., -20, 0.);                             // values are in kG
    fMagField->SetField(0., -4, 0.);                             // values are in kG
                                                                  //  (xmin,xmax,ymin,ymax,zmin,zmax)
     run->SetField(fMagField);
		// cout<<"By at (0,0,0) is "<<run->GetField()->GetBy(0,0,0)<<endl;
    // --------------------------------------------------------------------

    // -----   Create PrimaryGenerator   --------------------------------------
    FairPrimaryGenerator* primGen = new FairPrimaryGenerator();

    // Add a box generator also to the run
		 //some common PDGcodes:  2212(proton), 2112(neutron), 11(electron), 13(muon), 111(Pi0), 211(PiPlus), 1000010020(deutron) 1000010030(triton), 1000020030(He3),1000020040(He4),1000060120(C12)
		 FairBoxGenerator* boxGen = new FairBoxGenerator(2212, 3);   // 2212 = proton, 13 = muon; 1 = multipl.
    boxGen->SetPRange(0.4, 0.5);                                // GeV/c
    boxGen->SetPhiRange(0., 360.);                            // Azimuth angle range [degree]
    boxGen->SetThetaRange(0., 0.);                           // Polar angle in lab system range [degree]
		boxGen->SetXYZ(0., 0., -200.);                               // cm
		// boxGen->SetXYZ(0., 0., 0.);                               // cm
		
		// Int_t mult=2;
		// FairIonGenerator* boxGen = new FairIonGenerator(6, 12, 6, mult, 0., 0., 1.0, 0., 0., -200.);
		/***FairIonGenerator(z, a, q, mult, px, py, pz, vx, vy, vz);
		 **@param z         Atomic number
		 **@param a         Atomic mass
		 **@param q         Electric charge [e]
		 **@param mult      Number of ions per event
		 **@param px,py,pz  Momentum components [GeV] per nucleon!
		 **@param vx,vy,vz  Vertex coordinates [cm]
		 **/
    primGen->AddGenerator(boxGen);

    run->SetGenerator(primGen);
    // ------------------------------------------------------------------------

    //---Store the visualiztion info of the tracks, this make the output file very large!!
    //--- Use it only to display but not for production!
    run->SetStoreTraj(kTRUE);

    // -----   Initialize simulation run   ------------------------------------
    run->Init();
    // ------------------------------------------------------------------------

    // -----   Runtime database   ---------------------------------------------

    Bool_t kParameterMerged = kTRUE;
    FairParRootFileIo* parOut = new FairParRootFileIo(kParameterMerged);
    parOut->open(parFile.Data());
    rtdb->setOutput(parOut);
    rtdb->saveOutput();
    rtdb->print();
    // ------------------------------------------------------------------------

    // -----   Start run   ----------------------------------------------------
    run->Run(nEvents);

    // You can export your ROOT geometry ot a separate file
    run->CreateGeometryFile("geofile_full.root");
    // ------------------------------------------------------------------------

    // -----   Finish   -------------------------------------------------------
    timer.Stop();
    Double_t rtime = timer.RealTime();
    Double_t ctime = timer.CpuTime();
    cout << endl << endl;
    cout << "Macro finished succesfully." << endl;
    cout << "Output file is " << outFile << endl;
    cout << "Parameter file is " << parFile << endl;
    cout << "Real time " << rtime << " s, CPU time " << ctime << "s" << endl << endl;
    // ------------------------------------------------------------------------
}
